source("/home/jc227089/evo-dispersal/KBGrad/KBGradfunctions.R")
setwd("/home/jc227089/SRE/KB/DispOnly/limits/data")

fname<-"Hbyb_limits"
flist<-list.files(pattern=fname)
flist<-flist[!grepl(pattern="concat", flist)]

concat<-c()
poplist<-vector(mode="list", length=5*length(flist))
for (ii in 1:length(flist)) {
	load(flist[ii])
	sidelist<-vector(length=2, mode='list')
	for (jj in 1:length(out)){
		if (is.null(out[[jj]]$pop)) next
		poplist[[(ii-1)*5+(jj)]]<-out[[jj]]$pop
		xmin<--min(out[[jj]]$pop[,"X"])
		xmax<-max(out[[jj]]$pop[,"X"])
		pars<-unlist(out[[jj]]$parameters)
		Dbar<-mean(out[[jj]]$pop[,"DP"])
		Hbar<-mean(out[[jj]]$pop[,"HP"])
		out[[jj]]$ztrace<-na.omit(out[[jj]]$ztrace)
		for (ss in 1:2){
			rws<-which(out[[jj]]$ztrace[,'ss']==ss)
			mean.cline<-apply(out[[jj]]$ztrace[rws,1:6], 2, mean)
			var.cline<-apply(out[[jj]]$ztrace[rws,1:6], 2, var)
			#acf.cline<-apply(out[[jj]]$ztrace[rws,1:6], 2, function(x){sum(acf(x, plot=FALSE, lag.max=15)$acf>0.1)})
			concat<-rbind(concat, c(pop=(ii-1)*5+jj, rep=jj, side=ss, pars, xmin=xmin, xmax=xmax, Dbar=Dbar, Hbar=Hbar, 
				mean.cline=mean.cline, var.cline=var.cline))
		}
	}
}
save(concat, poplist, file="concat_Hbyb_limits.RData")


####################################################
### Figures ###

rs<-which(concat[,'side']==1 & concat[,'var.cline.asymm.l']<250000)
pdf(file="../figures/test.pdf")
		clrs<-as.numeric(as.factor(concat[rs, 'b']))
		plot(concat[rs,'xmin'], concat[rs,'mean.cline.dNx'], col=clrs)
		legend('topleft', legend=levels(as.factor(concat[rs, 'b'])), pch=21, col=1:max(clrs))
		plot(concat[rs,'h2D'], concat[rs,'mean.cline.dNx'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'mean.cline.wdx'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'mean.cline.asymm'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'mean.cline.asymm.l'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'mean.cline.dD'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'var.cline.dNx'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'var.cline.wdx'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'var.cline.asymm'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'var.cline.asymm.l'], col=clrs)
		plot(concat[rs,'h2D'], concat[rs,'var.cline.dD'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'acf.cline.dNx'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'acf.cline.wdx'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'acf.cline.asymm'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'acf.cline.asymm.l'], col=clrs)
		#plot(concat[rs,'h2D'], concat[rs,'acf.cline.dD'], col=clrs)	
		
		clrs<-as.numeric(as.factor(concat[rs, 'h2D']))
		plot(concat[rs,'xmin'], concat[rs,'mean.cline.dNx'], col=clrs)
		legend('topleft',legend=levels(as.factor(concat[rs, 'h2D'])), pch=21, col=1:max(clrs))
		plot(concat[rs,'b'], concat[rs,'mean.cline.dNx'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'mean.cline.wdx'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'mean.cline.asymm'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'mean.cline.asymm.l'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'mean.cline.dD'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'var.cline.dNx'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'var.cline.wdx'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'var.cline.asymm'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'var.cline.asymm.l'], col=clrs)
		plot(concat[rs,'b'], concat[rs,'var.cline.dD'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'acf.cline.dNx'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'acf.cline.wdx'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'acf.cline.asymm'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'acf.cline.asymm.l'], col=clrs)
		#plot(concat[rs,'b'], concat[rs,'acf.cline.dD'], col=clrs)		
		
dev.off()


limits<-concat[,"xmax"]+concat[,"xmin"]
zmat<-tapply(limits, list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
zmat2<-tapply(concat[,"mean.cline.dD"], list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
zmat3<-tapply(concat[,"mean.cline.asymm.l"], list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
zmat4<-tapply(concat[,"Dbar"], list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
zmat5<-tapply(concat[,"mean.cline.wdx"], list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
zmat6<-tapply(concat[,"mean.cline.dNx"], list(b=concat[,"b"], h2D=concat[,"h2D"]), mean)
#zmat[is.na(zmat)]<-0

pdf(file=paste("../figures/","ContourPlots", ".pdf", sep=""))
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat, xlab="b", ylab="h2D", main="Range limits: h^2_d by b")
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat2, xlab="b", ylab="h2D", main="Dispersal cline: h^2_d by b")
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat3, xlab="b", ylab="h2D", main="Asymmetry cline: h^2_d by b")
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat4, xlab="b", ylab="h2D", main="Mean dispersal value: h^2_d by b")
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat5, xlab="b", ylab="h2D", main="Mean fitness cline: h^2_d by b")
contour(x=as.numeric(levels(as.factor(concat[,"b"]))), y=as.numeric(levels(as.factor(concat[,"h2D"]))), 
	z=zmat6, xlab="b", ylab="h2D", main="Mean density cline: h^2_d by b")
dev.off()

# plot a pooled set of replicate populations for a particular value of b and H2D
b<-2
h2D<-0
i<-unique(concat[(concat[,"b"]==b & concat[,"h2D"]==h2D), "pop"])
 
pops<-poplist[i]
for (ii in 1: length(pops)){
	plotter.mean(pops[[ii]], K=150, lambda=15, b, stable.domain=50, bw=1, filename=paste("../figures/pop", ii, ".pdf", sep=""))
}

